home *** CD-ROM | disk | FTP | other *** search
/ Over 1,000 Windows 95 Programs / Over 1000 Windows 95 Programs (Microforum) (Disc 1).iso / 1249 / mat_util.t < prev    next >
Text File  |  1997-04-18  |  5KB  |  302 lines
  1. %
  2. % "mat_util.t" contains matrix utility functions for use
  3. % in other programs. 
  4. %
  5. % To use, add file to list in "*.prt" file
  6. %
  7. % Dimension 'DIM' should be declared as a global constant
  8. % in the main program file.
  9. %
  10. %   Sample program for the T Interpreter by:
  11. %
  12. %   Stephen R. Schmitt
  13. %   962 Depot Road
  14. %   Boxborough, MA 01719
  15. %
  16.  
  17. type rmatrix : array[DIM,DIM] of real
  18. type rvector : array[DIM] of real
  19. type ivector : array[DIM] of int
  20.  
  21. %
  22. % print a matrix
  23. %
  24. procedure print_mat( var m : rmatrix )
  25.  
  26.     var i, j : int
  27.  
  28.     for i := 0...DIM-1 do
  29.  
  30.         for j := 0...DIM-1 do
  31.  
  32.             put m[i,j]:15:6...
  33.  
  34.         end for
  35.  
  36.         put ""
  37.  
  38.     end for
  39.  
  40. end procedure
  41.  
  42. %
  43. % return the absolute value of a real number
  44. %
  45. function rabs( x : real ) : real
  46.  
  47.     if x < 0.0 then
  48.  
  49.         x := -x
  50.  
  51.     end if
  52.  
  53.     return x
  54.  
  55. end function
  56.  
  57. %
  58. % invert a square matrix using Gauss-Jordan method
  59. %
  60. % note: pivoting will sometimes fail; set pivot to false to disable
  61. %
  62. function invert( var a, b : rmatrix, pivot : boolean ) : real
  63.  
  64.     var w : array[DIM,2*DIM] of real
  65.     var p : array[DIM] of int
  66.     var s, t, det : real
  67.     var i, j, k : int
  68.     label invert_exit :
  69.  
  70.     % setup workspace : [A|I]
  71.  
  72.     for i := 0...DIM-1 do
  73.  
  74.         for j := 0...DIM-1 do
  75.  
  76.             w[i,j] := a[i,j]
  77.             w[i,j+DIM] := 0.0
  78.  
  79.         end for
  80.  
  81.         w[i,i+DIM] := 1.0
  82.  
  83.     end for
  84.  
  85.     det := 0.0
  86.  
  87.     % convert to upper triangular form
  88.  
  89.     for k := 0...DIM-2 do
  90.  
  91.         j := k
  92.         t := 0.0
  93.  
  94.         if pivot then           % find pivot row for column k
  95.         
  96.             for i := k...DIM-1 do
  97.  
  98.                 s := rabs( w[i,k] )
  99.             
  100.                 if s > t then
  101.  
  102.                     t := s
  103.                     j := i
  104.  
  105.                 end if
  106.  
  107.             end for
  108.  
  109.         end if
  110.  
  111.         p[k] := j               % save pivot row number
  112.  
  113.         if j ~= k then          % swap rows
  114.  
  115.             for i := k...2*DIM-1 do
  116.  
  117.                 t := w[j,i]
  118.                 w[j,i] := w[k,i]
  119.                 w[k,i] := t
  120.  
  121.             end for
  122.  
  123.         end if
  124.  
  125.     end for
  126.  
  127.     p[DIM-1] := DIM - 1
  128.  
  129.     % zeroize elements below each diagonal element
  130.  
  131.     for k := 0...DIM-2 do
  132.  
  133.         if w[k,k] = 0.0 then
  134.  
  135.             goto invert_exit
  136.  
  137.         end if
  138.  
  139.         for i := k+1...DIM-1 do
  140.  
  141.             t := w[i,k] / w[k,k]
  142.  
  143.             for j := k...2*DIM-1 do
  144.  
  145.                 w[i,j] := w[i,j] - t * w[k,j]
  146.  
  147.             end for
  148.  
  149.         end for
  150.  
  151.     end for
  152.  
  153.     % compute determinant
  154.  
  155.     det := 1.0
  156.     i := 1
  157.  
  158.     for k := 0...DIM-1 do
  159.  
  160.         if p[k] ~= k then
  161.  
  162.             i := -i
  163.  
  164.         end if
  165.  
  166.         det := det * w[k,k]
  167.  
  168.     end for
  169.  
  170.     det := i * det
  171.  
  172.     % make diagonal elements = 1
  173.  
  174.     for k := 0...DIM-1 do
  175.  
  176.         t := w[k,k]
  177.  
  178.         for j := k...2*DIM-1 do
  179.  
  180.             w[k,j] := w[k,j] / t
  181.  
  182.         end for
  183.  
  184.     end for
  185.  
  186.     % zeroize elements above each diagonal element
  187.  
  188.     for decreasing k := DIM-1...1 do
  189.  
  190.         for decreasing i := k-1...0 do
  191.  
  192.             t := w[i,k]
  193.  
  194.             for j := k...2*DIM-1 do
  195.  
  196.                 w[i,j] := w[i,j] - t * w[k,j]
  197.  
  198.             end for
  199.  
  200.         end for
  201.  
  202.     end for
  203.  
  204.     % output results
  205.  
  206.     for i := 0...DIM-1 do
  207.  
  208.         for j := 0...DIM-1 do
  209.  
  210.             b[i,j] := w[i,j+DIM]
  211.  
  212.         end for
  213.  
  214.     end for
  215.  
  216.  
  217.     invert_exit :
  218.  
  219.     return det
  220.  
  221. end function
  222.  
  223. %
  224. % multiply two square matrices
  225. %
  226. procedure mul_mat_mat( var a, b, c : rmatrix )
  227.  
  228.     var i, j, k : int
  229.     var s : real
  230.  
  231.     for i := 0...DIM-1 do
  232.  
  233.         for j := 0...DIM-1 do
  234.  
  235.             s := 0.0
  236.  
  237.             for k := 0...DIM-1 do
  238.  
  239.                 s := s + a[i,k] * b[k,j]
  240.  
  241.             end for
  242.  
  243.             c[i,j] := s
  244.  
  245.         end for
  246.  
  247.     end for
  248.  
  249. end procedure
  250.  
  251. %
  252. % multiply a square matrix by a vector
  253. %
  254. procedure mul_mat_vec( var a : rmatrix, var x, y : rvector )
  255.  
  256.     var i, j : int
  257.     var s : real
  258.  
  259.     for i := 0...DIM-1 do
  260.  
  261.         s := 0.0
  262.  
  263.         for j := 0...DIM-1 do
  264.  
  265.             s := s + a[i,j] * x[j]
  266.  
  267.         end for
  268.  
  269.         y[i] := s
  270.  
  271.     end for
  272.  
  273. end procedure
  274.  
  275. %
  276. % decompose a square matrix into symmetric and skew-symmetric parts
  277. %
  278. procedure sym_mat( var a, m, s : rmatrix )
  279.  
  280.     var i, j : int
  281.  
  282.     for i := 0...DIM-1 do
  283.  
  284.         for j := 0...DIM-1 do
  285.  
  286.             m[i,j] := ( a[i,j] + a[j,i] ) / 2
  287.  
  288.         end for
  289.  
  290.     end for
  291.  
  292.     for i := 0...DIM-1 do
  293.  
  294.         for j := 0...DIM-1 do
  295.  
  296.             s[i,j] := ( a[i,j] - a[j,i] ) / 2
  297.  
  298.         end for
  299.  
  300.     end for
  301.  
  302. end procedure